Critical decay exponent of the pair contact process with diffusion 
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We study the pair contact process with diffusion (PCPD) by extensive Monte Carlo simulations, 
focusing on the density decay exponent 5 at criticality. Studying how the ratio of the density and the 
pair density behaves at criticality, we estimate the exponent \ °f the leading correction to scaling as 
X = 0.37 ± 0.01. With \ so obtained, we study the effective exponent of the critical density decay, 
to conclude that S — 0.185 ± 0.010 which is different from the critical exponent of the directed 
percolation (DP), 0.1595. 

PACS numbers: 05.70.Ln,05.70. Jk,64.60.Ht 



I. INTRODUCTION 

The pair contact process with diffusion (PCPD) is an 
interacting particles system with diffusion, pair annihi- 
lation (2 A —> 0), and creation of particles by a pair 
(2 A — > 3 A) . A decade- long research has failed to elicit a 
full consensus concerning the universality class of the one 
dimensional PCPD (for a review of the early controversy, 
see Ref. [l[). Two competing hypotheses are struggling 
for survival; does the PCPD belong to the directed per- 
colation (DP) universality class or not 0-0]? 

The only consensus is that the PCPD has very strong 
corrections to scaling, which is the reason why reported 
numerical values of the critical exponents are widely scat- 
tered una. To overcome the strong corrections to 
scaling, a few approaches have been put forward. In 
Ref. [5], it was suggested that the so-called 'soft-core bo- 
son' PCPD model does not have such strong corrections 
and that the PCPD does not belong to the DP class. 
However, it turned out that the 'soft-core boson' model 
actually has strong corrections to scaling Q. Besides, 
the suggested exponent in Ref. 0] was disputed 0]. 

Rather than facing with the strong corrections to scal- 
ing, Park and Park suggested indirect tests about the 
universality class in a series of papers 0,0,111 El. When 
relative motion of particles and pairs is spatially biased, 
a completely different critical behavior from the DP as 
well as from the PCPD appears 0, [lij . Since this crit- 
ical change due to a spatial bias has not been observed 
in models known to belong to the DP class, this was 
argued to be the evidence that the PCPD does not be- 
long to the DP class 0|. When an interaction which 
clearly makes the system belong to the DP class is in- 
troduced to the PCPD, nontrivial crossover behaviors 
have emerged 0, [po| . Quite interestingly, a nontrivial 
crossover was observed between two models belonging to 
the DP class . On one hand, this nontrivial crossover 
in the DP class is mediated by the drastic shrink of the 
number of absorbing states (from infinity to one, so to 
speak). On the other hand, the crossover from the PCPD 



to the DP has nothing to do with the structural change 
of the absorbing space. Since the 'crossover' between DP 
models without structural change of absorbing space is 
mostly trivial (no exception has ever been reported), it 
is again claimed that the PCPD does not share critical 
behavior with the DP 0,[Hj]. 

Admittedly, however, the indirect tests cannot fully 
answer the question. Thus, we feel it necessary to study 
the PCPD extensively to find the critical exponent di- 
rectly. Although the strong corrections to scaling is a 
difficult obstacle hard to evade, it can be helpful if we 
know how strong the corrections are. In this paper, we 
carefully study the corrections to scaling. Then, using 
the information about the corrections to scaling, we will 
estimate the critical decay exponent. 

The paper is organized as follows: After introducing 
the model in Sec. UH we analyze simulation results in 
Sec. IIIII In Sec. IIV1 we discuss about our analysis and 
summarize our work. 



II. MODEL 

The model we are interested in is the one-dimensional 
parity-conserving pair contact process with diffusion first 
introduced in Ref. [171 ] . Since the parity conservation 
does not play any role in determining the universality 
class [lfl IT?} , we will simply refer to this model as the 
PCPD. The model is defined on a one-dimensional lattice 
of size L with periodic boundary conditions. The state 
at each site is either occupied by a particle (A) or empty 
(0); every site can accommodate at most one particle. 
The dynamics can be defined as 
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where a = (1— p) /2. The process starts from the fully oc- 
cupied initial condition. Obviously, a configuration with- 
out a particle is an absorbing state. Although a configu- 
ration with a single particle is not an absorbing state in 
a strict sense, we loosely define such a state as absorb- 
ing because a single particle alone in the system cannot 
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increase the number of particles. The PCPD in the con- 
ventional setting (2 A — > 0, 2A — » 3A) also has such two 
'absorbing' states. 

In simulations, we keep a list of the position of every 
particle for all time. Assume that there are N t particles 
at time t. Among Nt particles, we choose one with equal 
probability. Let us assume that the chosen particle is lo- 
cated at site i. With probability i, one of two nearest 
neighbors of site i is chosen. Let us call the index of the 
chosen site j (j is either i+ 1 or i — 1 ; recall that we are us- 
ing periodic boundary conditions) . If site j is empty, the 
particle at site i moves to site j. If site j is also occupied, 
two particles at both sites (i and j) are annihilated with 
probability p. But with probability 1 — p, the pair will 
try to branch another two particles. If a branching event 
is scheduled, we choose one of two directions with equal 
probability. If the direction is determined, we check if 
two nearest sites of the pair along the chosen direction 
are both empty. If so, two new particles are placed there, 
but otherwise nothing happens and there is no change of 
configuration. After the attempt, time increases by 1/N t 
and the list is appropriately updated. 

During simulations, we measured the (particle) density 
p(t) and the pair density p p (t) which are defined as 



P(t) 



~X>i(t)>, ft»(*) = ~X><(*><+i(t)), (2) 



where Si(t) is a random variable which takes 1 (0) if site 
i is occupied (vacant) and (...) means the average over 
ensembles. We also counted the number of runs which fall 
into one of the absorbing states, but all simulation results 
to be reported later did not end up with an absorbing 
state up to the observation time. 



III. CRITICAL DECAY EXPONENT 

At criticality, the density is expected to show power- 
law behavior with corrections to scaling such as 



p(t) ~ ai t- b (l + Clt~ X ) 



(3) 



where <5 is the critical decay exponent, x is the exponent 
of the leading correction, and a\, c\ are (nonuniversal) 
constants. For the DP class, 5 « 0.1595 in one dimension. 
Since the two dimensions are believed to be the upper 
critical dimensions of the PCPD, we assume that there is 
no logarithmic correction in the one dimensional system. 

To determine the critical exponent S from numerical 
data, it is customary to analyze the effective exponent 
defined as 



S e ff(t; b) 



\np{t) - \np{t/b) 
h76 : 



(4) 



where & is a time independent constant. At criticality the 
asymptotic behavior of the effective exponent should be 

-5 cS {t;b)^-8- Cl h ^—^t-^ + o{t-^). (5) 
in o 



Thus, if we plot the effective exponents — 5 e s as a func- 
tion of t~ x at criticality, the curve should approach to 
the critical exponent with finite slope as t~ x — > 0. On 
the other hand, if the system is in the active (absorb- 
ing) phase near the critical point, the effective exponent 
eventually veers up (down) as t~ x — > 0. The informa- 
tion about x, therefore, will be very useful to find the 
accurate value of the critical exponent. 

To find x without prior knowledge of S, we need an- 
other quantity which decays as t~ s in the asymptotic 
limit. A natural candidate is the pair density p p (t). Al- 
though a mean field theory and simulations of higher- 
dimensional model [HI as well as the one dimensional 
driven PCPD Q have shown that the pair density does 
not necessarily follow the same power law behavior as the 
density, numerical studies of the one-dimension PCPD 
(for example, see Ref. @, Q) show that p(t)/p p (t) indeed 
approaches to a finite number as t — ¥ oo at criticality. 
Hence, it is expected up to the leading correction that 



p P {t) 



a 2 t- & (1 



C 2 t 



(6) 



where a 2 and c 2 are constants. Assuming c\ ^ c 2 , the ra- 
tio of p(t) and p p {t) at the critical point can be expanded 
as 



R(t) 



P(t) 



P P {t) 



ax 
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[i + ( Cl - c 2 )r x (l + ct~ x )] , (7) 



where we have assumed that the next leading correction 
to scaling is not larger than t~ 2x and C ^ 0. This ratio 
has already been studied in Refs. 0, H| and, albeit in a 
different context, in Ref. [7[. Note that two assumptions 
have been employed to arrive at Eq. ([7]). The plausibility 
of these two assumptions will be discussed in Sec. IIVI 

Since S does not appear in R(t), we can study x directly 
using the effective exponent 



, LN , (R{t/b)-R{t/b 2 ) 



/In b, 



(8) 



where & is a time- independent constant. From Eq. (|7|), 
the asymptotic behavior of Xcs{t) at criticality is ex- 
pected to be 



Xef f(t; b)~ X + C 



b 2x - 1 
In b 



(9) 



Now we need the critical point. Actually, the critical 
point of the model in question is quite accurately esti- 
mated in Ref. [l5j . which was the reason why we study 
the present model rather than the conventional one. 

Equipped with Eq. @ together with the accurately es- 
timated critical point, we will find x by the Monte Carlo 
simulations. To this end, we simulated the system with 
L = 2 19 for several p's around the critical point (from 
0.180 210 to 0.180 215) up to t = 10 7 . Since there is 
no significant difference among these simulations up to 
t = 10 7 , we will only present the results for p = 0.180 214 
which later will be asserted to be the critical point. Be- 
cause we have to remove the leading behavior of R(t), 
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FIG. 1. (Color online) Plots of Xes{t; b) for b = 10 3/2 (upper 
curve), 10 2 (middle curve), and 10 5 ' 2 (lower curve) as a func- 
tion of t~ i7 . The straight lines are resulting fitting functions 
which were found as explained in the text. 



statistical noise could be enormous unless the number of 
ensembles is large enough. To reduce statistical noise, we 
collected data from 25 000 ensembles. Then, the effective 
exponents are calculated using b = 10 3 / 2 , 10 2 , and 10 5 / 2 . 

To extract \ from x g, we first analyzed Xeff(i; b = 
10 5 / 2 ) which is least noisy among three effective expo- 
nents. For a given x> we fit Xos(t', 10 5 / 2 ) using Eq. (|9]) 
with C to be a fitting parameter. We set the fitting range 
to be from t = 2 x 10 6 to 10 7 . Then, we graphically 
checked if Eq. ([9]) can fit the other effective exponents 
when the same C for given x is used. The best fit is at- 
tained when x — 0.37 is used; see Fig. [TJ In this case, C 
is estimated as ~ 2.15. Since x — 0-36 or 0.38 does not 
yield a nice fitting (see Supplemental Material at [URL] 
for the fitting results when x = 0.36 and 0.38 are used), 
we estimated x as 0.37(1) with the number in parenthe- 
ses to indicate the error. Note that this estimation of x 
is consistent with the previous studies [1, Hj] . In a sense, 
the leading correction to scaling is not terrifyingly huge; 
for t = 10 8 , the leading correction is of the order of 10~ 3 . 
Hence, with the present computing power, it seems fea- 
sible to find the critical exponent of the PCPD. 

Now we will move on to the analysis of the effective 
exponent <5 e ff using x — 0.37. We present the simula- 
tion results for p = 0.180 213 (300 runs), 0.180 214 (456 
runs), and 0.180 215 (300 runs) with larger system size 
[L = 2 21 ) and longer observation time (t = 10 9 ) than 
above. Up to t — 10 9 , no simulated system falls into 
an absorbing state, which minimally guarantees the neg- 
ligible finite size effect. Actually, all simulated samples 
always contain at least one pair. 

In Fig. [2J the effective exponents calculated with b = 
100 are plotted as a function of t~ x for different p's (see 
Supplemental Material at [URL] for the pair-density ef- 
fective exponents plot which looks more or less same as 
Fig.©. The data with p = 0.180 213 (0.180 215) veers up 
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FIG. 2. (Color online) Plots of — 5 c ft using b — 100 as a 
function of t~ x with X — 0.37 for p — 0.180 213 (upper curve), 
0.180 214 (middle curve), and 0.180 215 (lower curve). 

(down), which signals that the system is in the active (ab- 
sorbing) phase. In fact, the curvature for p = 0.180 215 
is not so conspicuous, but if we use b = 10, one can see a 
clear curvature for this case (see Supplemental Material 
at [URL] for the behavior of effective exponents calcu- 
lated with b = 10). Since the data with p = 0.180 214 
approaches to the ordinate with finite slope, we conclude 
that p c = 0.180 214(1) and 8 = 0.185(10). Our estima- 
tion of S should be compared with that of the DP class, 
0.1595. Quite interestingly, the critical decay exponent 
we found in this paper is close to the upper limit of S set 
by Hinrichsen Q. 



IV. SUMMARY AND DISCUSSION 

To summarize, we studied the pair contact process 
with diffusion with modulo 2 conservation. At first, we 
systematically measured the leading correction to scaling 
described by the exponent x- Then, by the analysis of 
the effective exponent 5 e g along with the estimated Xi 
we concluded that the critical density decay exponent is 
0.185(10), which is not compatible with that of the DP 
class. 

When we derive Eq. ([9]), however, we had to resort to 
a few assumptions. First, we assume that the coefficients 
of the leading corrections to scaling of the density and 
the pair density are different (c\ ^ c 2 ). Second (along 
with the first assumption), the next leading corrections 
to scaling is not larger than t~ 2x . We will discuss about 
these two assumptions in this section. 

The assumption of c\ ^ ci is actually not built on a 
firm ground. We cannot fully exclude the possibility that 
we actually estimated the exponent of the next leading 
correction, that is, c\ = c%. In fact, one can easily come 
up with two quantities whose coefficients of the leading 
correction are exactly the same. 

To see this, let us begin with writing the equation 



4 




where e is a (nonzero) constant. Hence if we plot 



FIG. 3. (Color online) Log-log plot of d(t) = (1 - p c )/(l - 
2p c )~ p p (t)/ p t (t) vs t at p = p c . A straight line segment whose 
slope is —1 is for the guide to the eyes. Inset: Log-linear plot 
of p p (t)/pt(t) vs t at p = p c . 



about the density behavior in time. A straightforward 
calculation gives 



-pit) = 2(1 - 2p)p p {t) - 2(1 -p)pt(t), 



where pt(t) is defined as 

t(*) = T "^(siSi+l [ s i+2 + (1 - S i+2 )s l+3 }). 



(10) 



(11) 



We drop the explicit t dependence of Sj for convenience 
and we have assumed that the system has mirror symme- 
try; that is, the density of the local configuration AA$A 
is assumed to be the same as that of A®AA. 

As in Sec. [TTT1 let us assume that at the critical point 
pit) and p p (t) behave in the asymptotic regime as Eqs. ([3]) 
and ([5]), respectively. Since Eq. ([TUf should be consistent 
with Eqs. and ([5]), Ptif) should behave up to the lead- 
ing correction as 



Ptif) 



a 3 t- s (1 



c 3 t 



(12) 



Since the leading behavior of the time derivative of p is 
, Eq. (fit)]) forces the relation (1 — 2p c )a 2 = (1— p c )d3- 



t 



-1-5 



By the same token, if \ < 1> c 2 should be equal to c 3 . 
In fact, any correction larger than t^ 1 should vanish on 
the right-hand side of Eq. (Till)) . Otherwise, the leading 
behavior of the term on the left hand side in Eq. (|TU|) 
which is t _1_<5 cannot be the same as the terms on the 
right hand side. Hence, if we study the ratio of p P (t) and 
Pt(t), it should approach to the asymptotic value as 



P P (t) _ 

Pt{t) l-2p 



l - Pc '1 + ei- 1 



oit- 1 )) 



(13) 



d(t) 



1 ~Pc 



Pp(t) 



(14) 



1 - 2 Pc p t (t) 
as a function of t at p — p c , one can see a 1/t behavior 
in the long time limit. In Fig. |3l we depict d(t) vs t in 
double-logarithmic scales at p — p c (while we were gath- 
ering data for the analysis of x> we also measured pt(t)). 
As a byproduct, clean 1/t behavior of d(t) supports that 
in Sec. IIIII we studied the regime where pit) and p p {t) are 
properly described by Eqs. (J3J) and ©. 

As the example of d{t) shows, it is not impossible that 
the leading correction to scaling in p and p p happens 
to be identical. Although we cannot provide a plausible 
argument, however, we think that Fig. [2] combined with 
Eq. ([5]) provides a numerical evidence that the leading 
correction of i?(t) is same as that of p. Note that similar 
studies in Refs. @, Q also implicitly assume that i?(i) 
indeed contain information about the leading correction 
to scaling. 

Actually, Fig. Q] which shows that Eq. ^ correctly 
describes the behavior of the effective exponent Xeff at 
the critical point also supports the assumption c\ ^ c-i. 
To see why, we will investigate what will be the cor- 
rect asymptotic behavior of XeB if c-x = c i- By this as- 
sumption, p and p p at criticality should take the form 
(Xi < X2 < 1) 

pit) = oi {1 + c<- Xl + ftt-^ + o(i-* 2 )} , 
P P (t) = a 2 {1 + ct-^ + f 2 t~™ + oit-^)} , (15) 
which gives 

02 R(t ^ 1 + Ct-K + fit-** 

ai ~ 1 + ci-** + f 2 t-*i 

«l + (/i-/ 2 )*- X2 (l-ct- XI )- (16) 

If this is the case, the procedure of finding \ explained 
in Sec. IIIII would not yield a reasonable conclusion. 

We can repeat the above disucssion to argue that the 
next correction to scaling is not stronger than t~ 2x . This 
assumption is also supported by Fig.Q]in a self-consistent 
way. Hence, two assumptions we made are supported 
self-consistently by the numerical study and the analysis 
in Sec. IIIII is legitimate. 
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SM FIG. 1. Fitting results of Xeff using X = 0.36 (left) and \ = 0.38 (right). 




SM FIG. 2. Effective exponent of the pair density with b = 100. 
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SM FIG. 3. Effective exponent of the particle density with b — 10. 
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